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A quasi-statistical equilibrium model is constructed to simulate the multicomponent composition 
of the crust of an accreting neutron star. The ashes of rp-process nucleosynthesis are driven by 
accretion through a series of electron captures, neutron emissions, and pycnonuclear fusions up 
to densities near the transition between the neutron star crust and core. A liquid droplet model 
which includes nuclear shell effects is used to provide nuclear masses far from stability. Reaction 
pathways are determined consistently with the nuclear mass model. The nuclear symmetry energy 
is an important uncertainty in the masses of the exotic nuclei in the inner crust and varying the 
symmetry energy changes the amount of deep crustal heating by as much as a factor of two. 
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I. INTRODUCTION 

A large set of observational data from accreting neu- 
tron stars, including X-ray bursts, superbursts, crust 
cooling, and the quiescent luminosity of transiently ac- 
creting sources has become available in recent years. This 
wealth of observational data leads to opportunities to 
probe the properties of neutron stars, reactions on exotic 
nuclei, and the nature of dense matter. However, the 
composition of the deepest layers of the crust of accret- 
ing neutron stars is still not yet well known. This work 
presents one of the first multicomponent crust models 
which computes the properties of the deepest regions of 
the crust. 

In an accreting neutron star, as accreted matter (prin- 
cipally hydrogen and helium) accumulates on the surface, 
nuclei in the outer crust are pushed to deeper layers. At 
densities near 10 6 g/cm 3 , the fusion of hydrogen and he- 
lium can become unstable and a thermonuclear explo- 
sion results, an X-ray burst. These X-ray bursts gen- 
erate heavier nuclei by burning hydrogen, referred to as 
rp-process nucleosynthesis pj. As matter accretes, the 
burst ashes are pushed to higher densities and undergo 
a series of nuclear reactions: electron captures, neutron 
emissions, and pycnonuclear fusions. These reactions 
drive the composition to nuclear statistical equilibrium, 
the ground state in the neutron star crust. This stable 
burning also generates heat of a few MeV per nucleon, 
which heats the crust in addition to the unstable burning 
which occurs in X-ray bursts. This is referred to as "deep 
crustal heating" Q . It is this deep crustal heating which 
is thought to drive the quiescent luminosity of accreting 
neutron stars 0|. At around 10 n ~ 12 g/cm 3 (the start of 
the inner crust) the neutron separation energy becomes 
negative and some of the neutrons form a quasi-free de- 
generate superfluid neutron gas. The neutron emissions 
dominate the deep crustal heating at densities just above 
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the density at which the inner crust begins. Finally, near 
10 14 g/cm 3 , the crust ends when nuclei are no longer 
energetically favorable. 

One success in connecting the observations to theo- 
retical models is the work in Refs. [1, [f|, which used a 
theoretical cooling model to describe the crust cooling of 
KS 1731 and MXB 1659. The X-ray flux from these ob- 
jects was observed immediately after outburst, and this 
flux decreased according to a broken power-law: the flux 
decreases more weakly with time at early times before 
the photons from the inner crust have reached the pho- 
tosphere, and the flux decreases more strongly with time 
at later times. This effect is due principally to the larger 
thermal conductivity from superfluid quasi-free neutrons 
in the inner crust (it is the first definitive observation 
of superfluidity in the crust). Ref. [j| showed that the 
crust must have a relatively small impurity parameter to 
have a thermal conductivity large enough to reproduce 
the data. 

Another success is the observation of the crustal cool- 
ing of SAX J1808 after an accretion outburst. SAX J1808 
is a neutron star with transiently accretes from a small 
main sequence companion. These accretion events warm 
up the neutron star crust relative to the core and the 
cooling of the crust can be observed right after the end 
of an accretion event. Ref. @ compared the observations 
to models developed in Ref. Q which required the com- 
position and heating of the accreted crust as an input. 
The result was that the cooling of SAX J1808 was so 
rapid as to suggest that the minimal neutron star cool- 
ing model H, Q was insufficient to explain the strong 
decrease in luminosity. This suggests that extra cooling 
beyond the minimal model, such as direct Urea or the 
cooling from pions or quarks, is present in SAX J1808. 

Some observations have been difficult to reproduce 
with current crust models. An instability to Carbon 
fusion at around 10 10 g/cm 3 may generate a super- 
burst [13, [ll|, an energetic form of X-ray burst which 
is observed in some sources. This instability is strongly 
temperature dependent. Current models suggest that 
the crust is too cold do destabilize Carbon fusion, thus 
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suggesting supcrburst models require revisiting [jj IT2I — 
[l4| . While the carbon fusion cross section is not well 
understood at these low energies, a severe enhancement 
in the fusion rate would be required to explain many su- 
perbursts [l5|. One important input parameter in these 
models is the amount of heating in the inner crust. If the 
heating in the inner crust was sufficiently strong, this 
could alleviate difficulties in current supcrburst models. 

Theoretical models of X-ray bursts which couple reac- 
tion networks to hydrodynamics and radiation transport 
have met with considerable success (l6| . These compli- 
cated computational frameworks cannot yet follow the 
reactions all the way the higher densities probed in super- 
bursts, in part because of the computational cost of the 
immense nuclear reaction networks which are required. 
An alternative was employed in Refs. [HI, [13], which uses 
a full reaction network with simplified hydrodynamics. 
These authors found extra heating from electron captures 
into excited states which partially, but not fully, allevi- 
ated the problem with supcrburst models. These sim- 
plified models become difficult in the inner crust, where 
nuclear reactions in the medium of quasi-free neutrons 
are not well known. Thus, these models cannot yet trace 
the evolution of accreted material as it sinks into the in- 
ner crust. A final difficulty is the requirement that the 
nuclear reactions are consistent with the nuclear masses, 
themselves subject to significant uncertainties at the rel- 
evant densities. 

On the other hand, multicomponent models of the ac- 
creted crust are particularly important because crusts are 
at a sufficiently low temperature and density to strongly 
limit the possible nuclear reactions which might occur. 
Nuclear structure effects also come into play, and some 
reactions which would not have been possible in a one- 
component system open up new reaction flows in a mul- 
ticomponent system. 

A simplification comes from the fact that, in the deep- 
est regions of the crust, most of the relevant nuclear re- 
action rates are strongly density dependent. This means 
that as soon as a reaction channel becomes energetically 
allowed, its rate rises and quickly becomes much faster 
than the local accretion timescale. A quasi-statistical 
equilibrium (QSE) ensues where, to a good approxima- 
tion, electron captures and neutron emissions always pro- 
ceed so long as they are exothermic. This QSE state is 
almost independent of the details of the nuclear reactions 
and depends more strongly on the masses of the nuclei 
which are present at any particular density (i.e. the Q 
values). Pycnonuclear fusion reactions are allowed when 
their fusion timescale is much faster than the local accre- 
tion timescale and can be handled separately. 

In this work, a one-zone multicomponent QSE model 
is used to describe the composition and heating in the 
crust of an accreting neutron star. When a particular 
nuclear reaction is energetically disfavored, its rate is pre- 
sumed to be zero. When non-fusion reactions are ener- 
getically favored, their rate is taken to be infinite, pro- 
ceeding until they become energetically disfavored again. 



Pycnonuclear fusion is handled by taking advantage of 
recent work on the relevant nuclear S-factors. Nuclear 
masses are described with a modern liquid droplet model 
which contains corrections for nuclear shell effects and 
matches available experimental data with an accuracy 
near that of more microscopic approaches. In-medium 
corrections are also included, i.e. the masses of each nu- 
clei depend on the temperature, the potential presence 
of quasi-free neutrons, and also on the ambient electron 
density. The use of a liquid droplet model allows the in- 
clusion of almost all of the relevant physics, but avoids 
the large computational time of more microscopic mod- 
els. 



II. CRUST MODEL 

This work assumes that the multicomponent crust is 
uniformly mixed, that nuclei are randomly distributed 
and uncorrelated, except for the lattice correlations 
which naturally occur in the Coulomb solid. The Wigner- 
Seitz approximation is used, and each nucleus occupies 
one and only one Wigner-Seitz cell, filled with electrons 
and quasi-free neutrons. Each cell is fixed in size by re- 
quiring it have no overall electric charge. 

The assumption of uniform mixing may be a poor one, 
especially at lower densities. Recent molecular dynam- 
ics calculations [l8l - [20j suggest that lighter and heavier 
nuclei separate. This will not affect the bulk energetics 
of the system but may have a significant impact on the 
pycnonuclear reaction rates described below. However, 
these microscopic simulations are computationally diffi- 
cult, and cannot yet be performed for a wide range of 
density regimes, nuclear mass models, and initial com- 
positions. These simulations also employ very simpli- 
fied models of the nucleon-nucleon interaction and do not 
handle all of the potential finite-size effects. 



A. In-medium nuclear mass formula 

The binding energy of the nucleus in medium with in- 
dex i can be written as a sum of terms, 

T) — -Ebulk + ^surf + 
^Coul + -Epair + ^shell ■(!) 

For simplicity, h = Hb = c = 1 in the following. The 
quantities n n .out and n Pj0 ut denote the average quasi- 
free neutron and proton densities outside of nuclei. The 
demarcation of nucleons inside and outside of nuclei is 
clearly artificial, but has proven to be a good approxi- 
mation for the bulk thermodynamic properties of matter 
except at the highest densities in the crust. The depen- 
dence of the mass formula on the n„ jC ,ut and n PiOU t, and on 
the ambient electron density, n e , is required to reflect the 
fact that the masses of nuclei depend on the surrounding 
medium. While the physical properties (radii, internal 
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neutron and proton densities, etc.) of each nucleus in 
the distribution are different for each nucleus and should 
thus be given a subscript i, this subscript is omitted in 
this section to simplify the notation. The eleven nuclear 
mass model parameters described below are all identi- 
cal for each nucleus, and the distinction between model 
parameters and nuclear properties will be made clear. 

The bulk part of the nuclear energy is constructed from 
a given equation of state (EOS) of bulk nuclear matter, 
denoted e^Un, n p ,T). Let £oo(n n , n p , T) denote the en- 
ergy density of bulk matter defined without the rest mass 
energy density, and similarly for the Helmholtz free en- 
ergy density, /, e.g. 

£ 00 (j^"n 1 Tip j -^J " £-00 {j^n 1 Tlpi -^ n ) T?!^/!^ TYlpTlp 
foo (j^n 1 Tip j T~) — foo {j^n ) ^p; TTl n Tl n THpTlp 

= Soo{n n ,n p ,T) -Ts(n n ,n p ,T) (2) 

The function e 00 (n„,n p ,T) is given either by a 
Skyrme [2l[ model or by the model of Akmal et al [22[ 
(hereafter APR). Skyrme models SLy4 [23j] . Gs, and 
Rs [24| are used, motivated by the fact that they provide 
a variation in the symmetry energy while still giving rea- 
sonable saturation properties and neutron star masses 
and radii. Finite temperature corrections in the bulk 
part of the nuclear energy are negligable at the tempera- 
tures of interest (< 10 9 K), so s(n n ,n p ,T) will be taken 
to be zero. In effect, the mass model described below 
is actually an entire class of mass models with different 
functions, Soo, all of which have a comparable quality 
as evaluated by the RMS deviation of the mass excess, 
yet with different compressibilities and symmetry ener- 
gies. This allows one to estimate the uncertainties in 
the property of the crust due to the uncertainties in the 
nature of the nucleon-nucleon interaction (25| . 

An approximate expression for the volume of the nu- 
cleus is Kuc = A/nB,ia, where n B . m = n nM + n pAn , the 
sum of the average neutron and proton densities inside 
the nucleus. Then the bulk energy is (c.f. [Ill), 
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wher 



and 



are the internal neutron and proton 



densities inside the nucleus. In the present model, the 
internal baryon density is chosen to be 



n B ,in = n + nil 



(4) 



where no and n\ are parameters of the model and I = 
1 — 2Z/A. The individual average neutron and proton 
number densities are given by 



n n> in = n B , in (l + 6)/2 + g(x) 
n p , in = n B ,m(l - 8)/2 - g(x) 



(5) 



and the density asymmetry S — £1 where £ is an addi- 
tional parameter of the model. Other models choose to 
treat n n ^ n and n p ^ n as nucleus shape parameters to be 



minimized over for each nucleus, but this does not typi- 
cally improve the fit to experimentally measured masses. 
If C = I , nuclei in vacuum have no neutron skin and £ < 1 
indicates the presence of a neutron skin. The fraction of 
the Wigner-Seitz cell volume occupied by the nucleus, Xi 
is described below. Without the additional correction, 
g(x)> the energy of ultra-neutron rich matter the densi- 
ties implied by n n ^ n and n Pi i n fail to give a physical value 
of the core-crust transition density as computed in, e.g. 
Ref. [1^]. The function g(x) is defined by 



g(x) = fnC 



1 



f 



(6) 



where f n E = 5 and f n c = 1/2 alleviates this diffi- 
culty. The nuclear radii are defined by the relations 
47ri?^n„, in = 3A and ^TrR^n p ; nl = 3Z. 

The radius of the Wigner-Seitz cell, i?ws for each nu- 
cleus is determined by assuming that each cell contains 
the same number of protons and electrons, i.e. 



ttR 3 



ws' 



Z 
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where the electron density n e is taken to be the same in 
the WS cells of all nuclear species. This choice approxi- 
mately ensures that the edges of every Wigner-Seitz cell 
is at a fixed electrostatic potential. Furthermore, this 
choice ensures that the WS cells for each nucleus in the 
multicomponent mixture occupy the entire volume and 
that each neutron in the quasi-free neutron gas is asso- 
ciated with one and only one WS cell. The volume frac- 
tion of the cell which is occupied by the neutrons in the 
nucleus is x = (Rn/Rws) 3 and the volume fraction oc- 
cupied by protons is Xv = (-Rj)/-Rws) 3 - Because the size 
of the cell depends on the ambient electron density, and 
because the Coulomb energy in each nucleus depends on 
the size of the cell, the mass of every nucleus depends (al- 
beit weakly) on the number density of every other species 
in the system. This is expected, since the Coulomb en- 
ergy has longer range than the nuclear forces and couples 
each nucleus to the others. Defining volume of each cell, 
Vi 

exactly holds 

The surface energy is given as 



4:TrRy VSi /3 ensures that the identity J2i n iVi = 1 



-Esurf = 47ri?f urf a'S(n„, n p ) 



4 3 
-7ri? surf n B 



where R SU rf is defined by the relation 

A, 

the quantity B is defined by 
B(n n ,n p 
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[l/x 3 + 6+l/(l-x) 3 ] 



(8) 
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and as = 96/(6+ 16). This is equivalent to Eq. 5 in 
Ref. [251 ] . and ensures that the surface energy properly 
obeys the a; 3 dependence shown in Ref. 27} to match 



Thomas-Fermi calculations of very neutron-rich nuclei. 
It has the consequence that nuclei at large densities have 
small surface energies because the vanishing proton frac- 
tion requires the neutron density distribution to be quite 
diffuse (even though the current model contains no ex- 
plicit diffusiveness). 
The Coulomb energy is 



Seoul = 2C7re R p (n pM - n PtOUt ) f d (x P ) 



where the function f c {Xp) is given by 
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where d is the dimensionality (shape) of the nucleus. All 
nuclear are assumed to be spherical, i.e. d = 3. The coef- 
ficient, C, is an arbitrary parameter which decreases the 
Coulomb energy slightly, in order to model the diffusive- 
ness of the proton distribution in laboratory nuclei. This 
formula can include pasta by allowing d to be different 
from 3, but this possibility is left to future work. This 
differs slightly from the original expression in p8l |: the 
factor (n Pi - m — n PjOUt ) ensures that the Coulomb energy 
vanishes if the proton densities internal and external to 
nuclei are equal. In practice, however, this correction will 
not affect our results. 
The pairing energy is 



-A-^ 3 E A N and Z even 
{ A-^ 3 E A N and Z odd 
otherwise 



(13) 



The exponent 1/3 is known to be not well constrained by 
fitting to laboratory nuclei, and varying this exponent as 
a model parameter does not substantially improve the fit. 

For the shell energy, the corrections described in 
Ref. are used, modified to correct for the medium. 
Shell corrections can be particularly difficult to evaluate 
in quantum mechanical models of nuclei at these densi- 
ties because of spurious shell effects which are generated 
by the boundary of the Wigncr-Seitz cell [30| ■ It is also 
unclear how to properly modify these shell effects for the 
medium. One expects that, as the number density of neu- 
trons outside of nuclei increases relative to the number 
density of neutrons inside nuclei, the shell effects become 
less pronounced. This is treated phenomenologically by 
applying quenching functions 



A, 



A, 



Tin. in ^-n,out 



(14) 



to the shell correction energy, i.e. 

Tl v TL v Zv^v . 

J2 = — ^ — A„ H — — A T , 



Quantity (Units) 


APR 


SLy4 


Gs 


Rs 


n (far 3 ) 


0.1786 


0.1789 


0.1479 


0.1504 


m (far 3 ) 


-0.1057 


-0.08760 


0.03355 


0.02920 


V 


0.8804 


0.8798 


0.8642 


0.8696 


a (MeV/fm 2 ) 


1.155 


1.154 


0.9772 


0.9906 


<T<5 


1.382 


1.251 


0.4446 


0.4597 


c 


0.8957 


0.8933 


0.9246 


0.9227 


E A (MeV) 


5.224 


5.226 


5.213 


5.218 


ai (MeV) 


-1.390 


-1.390 


-1.378 


-1.373 


a 2 (MeV) 


0.008931 


0.01001 


0.01300 


0.01264 


10 3 a 3 (MeV) 


2.380 


2.360 


1.865 


1.920 


a np (MeV) 


0.1133 


0.1137 


0.09897 


0.09944 


(Sjtirms (MeV) 


1.132 


1.124 


1.228 


1.200 



TABLE I. The model parameters and RMS deviation in the 
mass excess for the homogeneous equations of state used in 
this work. 
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A„A P 



(15) 



Then the final shell correction is 



E f 



hell 



aiS 2 + a 2 S% + a 3 S 3 + a np Snp (16) 



D z 



The neutron magic numbers are set to 2, 8, 14, 28, 50, 
82, 126, 184, 228, 308, and 406, as suggested by Ref. [p], 
and the proton magic numbers are all within the range 
accessible by experiment. 

In summary, the 11 free parameters in this model (out- 
side of the input equation of state of bulk nuclear matter 
which is a kind of parameter in itself) are the surface ten- 
sion in MeV/fm 2 , a, the surface symmetry energy 175, the 
correction factor to the Coulomb energy, C, the asymme- 
try parameter the central density parameters, no and 
ni which are expressed in units of fm~ 3 , the pairing en- 
ergy Ej\, and the four parameters for the shell effects, 
Oij 02, 03, et np . 

The results of the fit to the experimental mass data 
from Ref. [HJ are given in Table UlAl Note that, because 
of the inclusion of shell effects, the quality of the mass 
formula is about 1.2 MeV, much closer to the 0.7 MeV 
deviation observed for FRDM, and much improved from 
a typical liquid droplet model which has a deviation of 
2.6 MeV or more. It is also instructive to see how the 
parameters depend with the density dependence of the 
symmetry energy: APR and SLy4 have symmetry ener- 
gies which depend rather more weakly with density and 
Gs and Rs have symmetry energies which depend more 
strongly with density. It is clear that the surface symme- 
try energy is correlated with the symmetry energy, as 
is the parameter n\, but the pairing and shell parameters 
are only weakly correlated to the symmetry energy 
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B. Free energy of multi-component matter 

The free energy density of matter (without the rest 
mass energy density) is 

/ (i^i}' ^n,outj ^p,out; T) 

J2i [Enuc(Zi,Ni,n n ^ ut ,np^ u t,n e ,T)ni + fc(rii,T)] 
+ (1 - <p)foo(n ntOUt ,n P!OUt ,T) + f c \ cc {n ei T) (17) 

where fc is the classical expression for the free energy 
density of nucleus i with density rii at temperature T and 
/dec is the free energy density of electrons without the 
electron rest mass energy density. One can also in princi- 
ple include finite temperature corrections to fc appropri- 
ate to dense matter as described in Ref. [33| . The contri- 
bution from fc is small and will be omitted here. Other 
thermodynamic quantities, such as the Gibbs energy den- 
sity (see the Appendix), can be trivially obtained from 
the free energy density in the usual way. Hereafter, the 
proton drip is assumed to be negligable, i.e. n PjOU t = 0. 
The electron density is not independent, n e must be self- 
consistently determined from the nuclear densities, i.e. 
n e = riiZ{. The partial volume fraction 

i i 

is also not independent and is a function of {rii}. In 
this formulation, n ritOU t is the local number density of the 
quasi- free neutron gas in between nuclei, while n.n,out(l — 
</>) is the number of neutron per unit volume inside a large 
volume of many Wigncr-Seitz cells. In a one-component 
system, tj) — \- 

The rest mass part of the energy density (not included 
above) is 

P = Ercst = (Nim n ni + ZiUiprii) + 

i 

(1 - 4>) n n ^ out m n + m e n e . (19) 

Note that the rest mass energy density is defined in terms 
of neutron and proton degress of freedom, even though 
the protons are typically bound in nuclei, which is con- 
venient since nuclear binding energies are being modified 
by the medium. 

The excluded volume correction can be written 

fexe = <Pifoo(n n ,out, 0, T) (20) 

i 

One can think about this either as a correction to the 
total energy density as is written in Eq. [T7] or as a cor- 
rection to the bulk energy of each nucleus 

^exc.i = - f — ) /oo(n-n,out,0,T) . (21) 

When written this way, one is effectively defining the 
mass of the nucleus in the medium relative to energy of 



pure neutron matter in an WS cell of equivalent volume, 
rather than defining the nuclear mass relative to the vac- 
uum. (This is similar to what has been historically done 
in some Hartree-Fock and Thomas-Fermi calculations in 
order to remove spurious shell effects [HI, HH) If one de- 
fines 

N ! = N . _ ^ Ihh^ =Ni _ 4 Rl n (22 ) 

3 

and define £ n uc,j = E nuc i + i? C xc,i then the free energy 
density of matter can be rewritten 

/({";}, »Vout,0,T) = 

(^nuc,i 
H~ foo (^n,out ? 

0,T) + / clcc (n e ,T). (23) 

This form now explicitly includes the rest mass part of 
the energy density and also makes computing analytical 
derivatives of the free energy a bit simpler. Note that 
some authors also refer to finite-volume corrections to 
fc(ni,T) as excluded volume corrections, but these cor- 
rections are negligable here. The total baryon density 
is 

n B = /^2 A%n ' 1 + ( 1 _ 4>) n n,out (24) 

i 

This formulation (with some minor additional finite- 
temperature effects) of the properties of dense matter is 
rich enough to express the properties of stellar matter 
in thermodynamic equilibrium at higher densities, con- 
ditions relevant for Type II supernovae. In this case, the 
multicomponent nature has a small impact on the over- 
all thermodynamic properties (27j thus justifying a "sin- 
gle nucleus approximation" pioneered in Ref. [36[ where 
matter was assumed to consist only of neutrons, protons, 
alpha particles, and representative heavy nuclei. Multi- 
component calculations of an EOS for supernova simu- 
lations have been performed in several works [H, l37l44~fl | 
but these works have not addresed the matter in the ac- 
creted neutron star crust. 



C. Quasi-Statistical Equilibrium 

In lieu of a full reaction network, nuclear reactions pro- 
ceed in "chunks" ; a small chunk of the nuclei present in 
the current distribution arc assumed to instantaneously 
undergo a particular reaction at constant pressure and 
constant entropy. In the case that the Hclmholtz free 
energy (enthalpy) per particle is lowered, then the re- 
action proceeds and the composition is modified accord- 
ingly. At the temperatures of interest, finite-temperature 
effects are negligible so the enthalpy is replaced with the 
Gibbs energy. Minimizing the Gibbs energy per particle 
at fixed pressure rather than the free energy density has 
the additional benefit that prevents the system from us- 
ing nuclear reactions to unphysically lower the pressure 
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as the baryon density is increased. As the chunk size ap- 
proaches zero, this procedure is equivalent to choosing a 
quasi-cquilibrium state at each density. The chunk size 
is chosen to be l/100th of the total number density of all 
nuclei, and this size is sufficiently small to ensure that 
the results approximate the correct quasi-equilibrium. 
This process does not a represent true statistical equilib- 
rium because the pycnonuclcar fusion reactions are not 
reversible, i.e. fission is not allowed. 

In order to match an old and new configuration at 
constant pressure, the nuclear densities and the global 
quasi- free neutron density, n„ iOUt (l — </>) in the new config- 
uration are scaled by the same factor (denoted a below) 
until the pressure matches that of the original configu- 
ration. Neutron emissions and neutron captures require 
an additional constraint to ensure that the neutrons are 
counted correctly. Denoting Sn as the decrease in the 
number density of parent nuclei due to neutron emission, 
the number density of parent nuclei after neutron emis- 
sion is a(np — Sn) and the number density of daughter 
nuclei after neutron emission is a(riD + Sn). Denoting the 
new configuration with apostrophes, The number density 
of quasi-free neutrons is given by 



t (l -(j>') = a [n„,out(l -<j>) + 6r 



(25) 



where a is the volume scaling factor and <j>' is implicitly 
a function of n' n out . The variable n' n out must be varied 
to ensure that this equation holds. In combination with 
the requirement that the pressures are equal, this gives 
two equations which must be solved for the variables a 
and n' n out to compute the proper new configuration af- 
ter any proposed nuclear reaction. This procedure ap- 
plies the neutron emission before the compression, but is 
equivalent to the opposite choice, i.e. n' P = anp — Sn, in 
the limit that Sn is small. 

The simulation of a crust begins with an initial com- 
position at a density near 10 6 g/cm 3 and proceeds in a 
series of quasi-equilibrium configurations to higher den- 
sities. Electron captures, /3-decays, neutron emissions, 
and neutron captures are always allowed if they lower 
the Gibbs energy per particle. In practice, since matter 
is becoming more neutron-rich as the density increases, 
electron captures and neutron emission dominate over j3- 
decays and neutron captures. Nuclei with number densi- 
ties less than 10 10 times the total number density of nu- 
clei are automatically pruned from the distribution. The 
heating rate is determined by computing the change in 
Gibbs free energy per baryon as the system proceeds from 
one configuration to another. In order to estimate the 
reduction in heating from neutrino emission in electron 
captures, heating from electron captures is multiplied by 
1/4 as in Ref. j42|. Because the liquid droplet model is 
not very accurate for very light nuclei, all reactions which 
result in nuclei with Z < 4 are not permitted, but this 
is a good approximation throughout the crust, as shown 
below. Nuclei whose neutron or proton radii are larger 
than the size of the WS cell are unphysical and thus also 
disallowed. 



D. Pycnonuclear Fusion Reaction Rates 

Pycnonuclcar fusion reaction rates are the most uncer- 
tain rates involved in the accreted crust, in particular the 
effect of large neutron skins on rates is not clear . Pre- 
viously, detailed fusion rates have not been widely avail- 
able, for example, in 0], a simplified fusion rates were 
based on approximate S factors obtained from the pa- 
rameterization in [43j . This situation has been improved 
two-fold: (i) fusion S-factors for light neutron-rich nu- 
clei were computed in [44| . and (ii) a detailed formalism 
for fusion rates in a multicomponent plasma has been 
described in [45| . 

Ref. 44[ computes S-factors for Z = 6, 8, 10, and 12 
isotopes with even neutron numbers. In the work below, 
it is assumed that Z=4 nuclei always fuse, though this 
assumption does not significantly affect our calculations. 
Nuclei Z=14 are assumed to fuse whenever Z=12 nuclei 
fuse, and this assumption is tested below. In order to 
compute S-factors involving odd proton or neutron num- 
bers, the proton number is always increased by one and 
the neutron number is decreased by one, which slightly 



decreases the fusion rates. Ref. [44| does not compute 
S-factors for nuclei which are sufficiently neutron-rich 
for our study, so it is assumed that the S-factors for all 
Z = 6 nuclei with A > A max (Z) are the same as that for 
A = A max (Z), where A max {Z) = 24, 28, 40, and 46, for 
Z = 6, 8, 10, and 12, respectively. This choice slightly 
increases the potential for fusion. 

To compute the rates, the formalism in Ref. [45| is 
used, which assumes that the multicomponent plasma 
is uniformly mixed. To compute the reaction rates, the 
parameter A is defined with 



As + A, 



AiAjZiZj [Z- /3 + zy s 



p(l-X n )(Z) 



(A) 1.3574 x 10 n g cm" 



1/3 



and the plasma temperature 



f 4irZiZje 2 nij\ 1 



p V 2 My 



where 



7T [(47m,/3r 1/3 + (47rn i /3) _1/3 



(26) 



(27) 



(28) 



for nuclei i and j with number densities rij and rij, and 
the reduced mass = m u AiAj / (Ai + Aj). With these 
definitions the rate is given by 

8p(l-X n )x i x j A i A j (A)ZfZj 
(1 + S^iAi + Ajf 



x S(E pk )\ 



3— C p i 



exp 



C e 



cm^s" 1 (29) 
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FIG. 1. The composition of the equilibrium neutron star crust 
as a function of density. Deviations originate primarily be- 
cause of the symmetry energy: models where the symmetry 
energy depends more steeply with density have larger Z and 
smaller N, i.e. a composition closer to the valley of stability. 
The smoother curves labeled "no shell" do not include shell 
elfects. 



In order to test if a particular fusion reaction is allowed, 
the the fusion timescale ni/R pyc is compared with the 
accretion timescale y/M. If the accretion rate timescale 
is larger, than the pynonuclear reaction is allowed if the 
fusion will lower the free energy, otherwise, the fusion is 
prohibited. 



III. RESULTS 

It is useful to compare the accreted crust to the much 
simpler case of an isolated neutron star crust in full equi- 
librium (all possible nuclear reactions are allowed) at zero 
temperature. The composition in the equilibrium crust 
is given in Fig. [TJ The neutron magic numbers are in- 
dicated by dashed lines, and it is clear that shell effects 
dominate the choice of the neutron number over all densi- 
ties, while the proton number steadily decreases to select 
the proper equilibrium electron fraction as it decreases 
with density. The proton number is most often even, be- 
cause of pairing. Note that at these small temperatures, 
the single nucleus approximation is very good because 
the energy difference between nuclear masses is typically 
much larger than the temperature and because the sys- 
tem is in equilibrium. 

The free energy per baryon of matter in the full equi- 
librium is given in Figure [2] Note that the uncertainty in 
the homogenous matter EOS affects the free energy per 



baryon at higher densites even though the composition 
is not strongly affected. 

Fig. [3] summarizes the properties of the accreted crust 
for the four EOSs used in this work. The entire crust is 
fixed at a temperature of T=10 8 K, and thermal effects 
do not strongly affect the composition. The initial com- 
position is taken from Table I in Ref. 18[. The lower 



panel for each model shows the average neutron and pro- 
ton numbers in the crust as a function of density along 
with the impurity parameter, (Q). Generally, neutron 
numbers increase with density and proton numbers de- 
crease with density, but a large amount of variation is 
present in the deepest regions in the crust. The dotted 
line in the upper panel is the chemical potential of the 
quasi-free neutrons in homogeneous matter at the same 
number density. The dashed line is the baryon chemical 
potential of the full equilibrium crust. The solid line is 
the full baryon chemical potential of matter in the ac- 
creted crust. 

The dashed-dotted line is the sum of this baryon chem- 
ical potential with the total integrated heat from all the 
layers above the current one. This quantity demon- 
strated in Fig. 5 of Ref. [42| to help explain the deep 
crustal heating phenomenon. The total amount of deep 
crustal heating is limited to the energy difference between 
this value and the baryon chemical potential in the full 
equilibrium crust. In models SLy4 and APR, the amount 
of deep crustal heat is limited by the rise of the baryon 
chemical potential in the full equilibrium crust. In mod- 
els Gs and Rs, where the symmetry energy depends less 
strongly with density, the baryon chemical potential of 
the full equilibrium crust is smaller and thus the deep 
crustal heating is more significant. 

Variations in the initial composition of matter are ex- 
plored in the top part of Fig. 2] The left panel has an 
initial composition of pure 106 Pd and the right panel has 
an initial composition of pure 56 Ni. The final compo- 
sitions in these models are not strikingly different from 
the X-ray burst ashes used in Fig. [3] The impurity pa- 
rameters are much lower, except for a peak near 4 x 10 12 
g/cm 3 in the 106 Pd case where the composition is re- 
covering from the larger number of neutrons stored in 
nuclei in that case. Increasing the temperature to 10 9 K, 
or varying the mass accretion rate (important for com- 
puting the relevant accretion timescale to compare with 
the fusion timescale) by an order of magnitude did not 
significantly change the results. 

The majority of the heating occurs through the emis- 
sion of neutrons near neutron drip. This heating often 
occurs in large chains which effectively convert some nu- 
clei entirely into neutrons, leaving the relative composi- 
tion of the remaining nuclei unaltered. A sample chain of 
this form begins with two 40 Mg nuclei, both of which un- 
dergo 6 electron captures and enough neutron emissions 
to form 22 C. These two nuclei fuse to form one 44 Mg nu- 
cleus, which then may emit 4 neutrons to return to the 
original 40 Mg. This cycle effectively converts one nucleus 
entirely into neutrons. Only a fraction of nuclei can un- 
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FIG. 2. The energy per baryon of matter in the cold neutron 
star crust in full equilibrium. The distinction between the 
APR and SLy4 models and the Gs and Rs models is primarily 
due to the difference in the density dependence of the nuclear 
symmetry energy. 



dergo such a cycle at any depth, and thus such cycles are 
difficult to resolve in single-nucleus models. 



IV. DISCUSSION 

This work is an important tool to calibrate more so- 
phisticated network calculations which follow the evolu- 
tion of nuclei in the crust as they begin at the surface 
and evolve to lower depths. 

Compositional nonuniformity may also drive novel 
heating processes in the crust 0, |4?| ■ While these heat- 
ing processes are not included in the current model, they 



are also fundamentally limited by the difference in the 
baryon chemical potential between the accreted and equi- 
librium crusts. 

The model outlined in this work ignores proton and 
light- fragment emission in accreting neutron stars, and 
may be incorrect near the crust-core phase transition. 
Near this density, protons may tunnel between nuclei and 
this will also affect the composition. Finally, the presence 
of pasta structures and nuclear structure uncertainties in 
the masses of nuclei are important in the deepest layers 
of the crust. The amount of heat generated at this depth 
is small, and the composition at the deepest regions in an 
accreting neutron star may be difficult to observe. The 
composition of the cold crust in equilibrium might be 
observable in the giant flares emitted by magnetars [48| . 

A one-zone model may not properly estimate the prop- 
erties of matter just near neutron drip. To demonstrate 
this, consider Fig. 5 of |42j]. In the single- nucleus ap- 
proximation, the baryon chemical potential, drops 
discontinuously as a function of the pressure at loca- 
tions where the composition of the crust changes. This 
is clearly unphysical, as this means that the system can 
lower its energy by a finite amount by taking a neutron 
from pressures slightly lower than the drop in /ib to pres- 
sures slightly larger than the drop. This means that neu- 
trons might sink as they are emitted, and the crust will 
reconfigure itself to ensure that there are no strong dis- 
continuities in fiB- (The electron chemical potential may 
also exhibit discontinuities in regions where multiple elec- 
tron captures occur.) The multicomponent models in this 
work soften these discontinuities considerably, but a flow 
of neutrons cannot be ruled out in this one-zone model. 
A multi-zone generalization of this work is in progress. 
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VI. APPENDIX: GIBBS FREE ENERGY 
DENSITY 

To compute the Gibbs free energy density, it is useful 
to compute the chemical potentials and entropies analyt- 
ically from Eq. |T7J An advantage of the liquid droplet 
model formalism is that these expressions are simple to 
compute accurately. 

The global quasi-free neutron density is defined with 
n n = n n ,out(l — 0), and this (not n„. ut is the quantity 
directly connected to baryon number conservation. Thus, 
to compute the Gibbs energy the derivatives 

K = (if) (30) 
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and 



df_ 

dn. 



(31) 



1 / n n ,{n,j'ij=£i},T 



are required. The symbol v is used rather than /i which 
is reserved for the chemical potentials in infinite matter 



/oc("n,out,T) 



/clcc (n e ,T) 



and 



(32) 



Note that these chemical potentials are defined including 
the rest mass energy. To be more concise, the subscripts 
will sometimes be suppressed in the following. The Gibbs 
free energy density is then 



(33) 



It is useful to note that 7i n ,m,i = ?in,m,i(Xi) is distinct 
quantity for each nuclear species i and is a function of 
Xi- For simplicity, this is denoted as n n i below. Also 
note that no separate term in the Gibbs free energy is 
required for the electrons since their number density is 
not an independent variable. In this section, it is also 
easier to think of the function £ nuc as a set of functions 
of the form £ n uc,i(n n ,out, Xi) ■ For each nucleus i, 



/ , n i Z 3 



ZiXi n r, 



(34) 



and this implies that \ an d n e both depend only on the 
number densities of nuclei and not separately on n„ iOUt . 
From Eq. 1341 one can obtain the relation dn e / ' drii = Zi 
and the derivative 



which shows that 0j is constant when all of the {rii} are 
held fixed. Thus, 



dn, 



dh r , 



{»«} 



and 



dn n . 



dn; 



1 ' {"jW/i, 



i - cb? 4" 



dm 



(36) 



(37) 



One can rewrite the derivative of <j)j in terms of the 
derivative of \j determined above 



9<Pj 
dn, 



N; 



' ii.l 



_ "jff'(Xj) ( dXj 
dn; 



' »j 



(38) 

where 5y is 1 if i = j and zero otherwise. The last 
derivative which is required is 



dn; 



-A-kBI 



+n r 





'dn n , out ^ 


3 ( 


y dn, j 


dR n j 




dXj 





(39) 



The derivative of / with respect to n n is 



dn n . 



dn„ 
An 



fd£ a 



\dn r , 



Note that this implies 

nn V Un ^n.out 



df 



(40) 



(41) 



The derivative of / with respect to the number density 
of nucleus i is 

df\ ( dn n , out \ 

Vi = [ I I I + ZiH e + fcnuc.i + Zim p 



<9n„ iOUt / \ dni 
+N-m n + njhi 



(42) 



dXj 
dn; 



N 3 Z t 



{rij}Vj/i, n„ 



Zjn nj 



z iXj9'{Xj) 



(35) 



This derivative is non-zero for i ^ j because the nuclei 
are all coupled by the Coulomb interaction through the 
electron density. Some simplification is afforded by the 
relation tf>i = Nini/n n i, where n n i depends only on 



where hij is defined by 
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